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Ecological systems comprise an astonishing diversity of species that cooperate or compete with 
each other forming complex mutual dependencies. The minimum requirements to maintain a large 
species diversity on long time scales are in general unknown. Using lichen communities as an 
example, we propose a model for the evolution of mutually excluding organisms that compete for 
space. We suggest that chain-like or cychc invasions involving three or more species open for creation 
of spatially separated sub-populations that subsequently can lead to increased diversity. In contrast 
to its non-spatial counterpart, our model predicts robust co-existence of a large number of species, 
in accordance with observations on lichen growth. It is demonstrated that large species diversity 
can be obtained on evolutionary timescales, provided that interactions between species have spatial 
constraints. In particular, a phase transition to a sustainable state of high diversity is identified. 



Introduction - Interactions between biological species 
may well be as old as life itself [IH3j with competition 
and predation as major determinants for species diversity 
[H [S] . Competitive exclusion [SI [7] has been suggested to 
reduce ecosystem diversity when several species compete 
for the same resources. Real ecosystems, on the other 
hand, consist of multiple species and have a robustness 
that may even increase with diversity [3 H] • To obtain 
robustness of ecosystem diversity in theoretical models 
one needs first of all to limit the exponential growth by 
assuming a maximum carrying capacity for the popula- 
tion of each species [10]. Extreme version of such mod- 
els [m [Hj indeed predicts a sustainable but fragile co- 
existence of multiple species. 

A more robust way to maintain high diversity is to in- 
clude space [T5HT5] e.g. in combination with hypercycles 
[m or predator-prey cycles 20-23 . Studying ecosys- 
tems in marine hard-substrate environments Jackson & 
Buss 24 suggested that non-transitive allelopathic rela- 
tionships between species could maintain species diver- 
sity on a longer timescale than pure hierarchical preda- 
tion relationships. This was confirmed in a model on 
two-dimensional lattice where sessile species compete for 
space [25l |26]. As in the "Buss" model [25l [26] we con- 
sider a 2-dimensional lattice where the competition for 
resources is a zero sum game about available space. In 
our model, however, the focus is on the dynamic balance 
between an introduction of new species and exclusion of 
older species. With this complementary model we, for 
the first time show that: a) There is a sharp transition 
from multiple to single species as the number of inter- 
actions is increased, b) Both cycles and chain-like (hi- 
erarchical) relationships lead to spatial fragmentation of 
species population thus creating isolated niches for new 
species and increased diversity. 

The model is inspired by the spatial dynamics of 
lichen communities. Lichens have existed for as long 
as ~600 million years [22 and are organisms consisting 
of fungi and algae living in a symbiosis |28j . Commu- 
nities of lichens are formed by a combination of slow 




FIG. 1. Photograph of a crustose hchen community on a rock 
in an alpine environment (at 1300m altitude, Jotunheimen, 
Norway). 



local growth and a reproductive strategy where fungal 
spores or propagules containing the intact symbiosis are 
dispersed (e.g. by winds or water flows) over lengthscales 
much larger than the size of the communities [29]. Be- 
cause of this long distance dispersal, we assume that new 
species come from far away, and are completely unrelated 
to any of the species that already colonize our system. 

Fig. [T] shows a crustose lichen community in an alpine 
environment. These lichens grow about O.lmm/year and 
typically covers a rock surface that has recently been ex- 
posed on a timescale ~ 100 — ?► 1000 years [50] • When 
a crustose lichen meets another, a contact boundary is 
formed, and if they are competitively equal the bound- 
ary remains stable over time. The bulging boundaries 
between various species seen in Fig. [T] suggest a mini- 
ecosystem with complex interactions. The interaction of 
these species may be represented by a directed network 
with directions from superior to inferior species. This 
network does not necessarily have any particular species 
as the most fit one in agreement with the fact that any 
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FIG. 2. Visualization of six successive snapshots of a simula- 
tion with 7 = 0.1 and a = 0.1 for a system of size L = 200. 
The networks below the snapshots illustrate the actual (green ) 
and potential (grey) links between the species in the system. 
The size of a node represents the current population size of 
the corresponding species whereas the thickness of green links 
quantifies the number of active invasions sites. 

sizable rock typically hosts multiple lichen species. 

Results. - Our model considers a multi-species ecosys- 
tem competing on a two dimensional lattice of sites which 
at any given time can be occupied by one species only. 
The species on one site can invade a neighbor site, pro- 
vided that it is occupied by a competitively inferior 
species. The emerging ecosystem can be characterized 
by a directed network of possible species interactions. 
Interactions which are materialized only when organisms 
of the respective species are neighbors somewhere in the 
system. The aim of our model is to study ecosystem di- 
versity as we change the number of potential interactions 
between species, parametrized by 7. In addition to this, 
new species are introduced with a rate a. Each time step 
of our model consists of two possible events: (i) Select a 
random site i and one of its nearest neighbors j. If the 
species s(i) at site i can invade the species s(j) at site 
j, i.e. r(s(i), s(j)) = 1, then site j is updated by setting 
s(i) = s(i). Here F is the matrix that represents the pos- 
sible interactions. These interactions remain fixed once 
they are introduced, (i) With probability a x 7/7V a new 
random species s is introduced at a random point j and 
assigned random interactions F(s, u) and T(u, s) with all 
existing species u in the system. Each of these interac- 
tions are assigned value 1 with probability 7, or otherwise 
set to (we do allow for the case T{u, s) = F(s,m) = 1). 
The introduced species s is assumed to be able to invade 




FIG. 3. Active interactions between species in the system. 
The number of interactions are measured in units of the max- 
imum potential interactions (sum of all outgoing links in the 
network) and shown as function of 7 and a. The computa- 
tions were carried out on a system of size L = 200. 



the previous species at the site j, s{j): F(s,s(j)) = 1. 

In Fig. [2] we show snapshots of a model ecosystem of 
size N = L X L = 200 x 200 using open boundary con- 
ditions. The snapshots represent a typical steady state 
behavior using 7 = 0.1 and a = 0.1. The different col- 
ors represent different species. Parts of the lattice are 
inactive (e.g. the orange colored species in the upper 
right corner), while other parts are exposed to active cy- 
cles, like the 4 mutually invading species grey — >■ dark 
brown — > brown — > red — > grey. This cycle is also illus- 
trated in terms of a network in the lower panels where the 
thickness of green links represents the number of active 
invasion events per time unit. 

With time the diversity (number of species in the sys- 
tem) D reaches a steady state where there is a balance 
between new species eliminating more than one of its prey 
species, and an increase in D when a new species only 
invades a fraction of an existing species. An increase in 
diversity requires that the population of existing species 
are fragmented into spatially separate regions. Two lo- 
cally interacting species cannot increase spatial hetero- 
geneity: One species will eventually be replaced by the 
other. To generate spatially separated regions, one needs 
simultaneous dynamics of at least 3 species locally. Such 
multiple interactions may be both chain-like or cyclic. 
For example a species A that emerges inside the territory 
of B, can result in a fragmentation of B if a third species 
C can invade A but not B. Whether B invades C or not, 
distinguishes a cycle from a chain-like relationship, but 
both can generate spatial fragmentation. However the 
cyclic relationship tends to live longer. 

In the network F a link is called active if it represents 
an interaction between two species which are physically 
in contact. For low 7 only a small fraction of the links 
F(s,m) = 1 are active since the species s and u are of- 
ten physically separated by other species which neither 
can invade. For larger 7 we observe fewer and more 
widespread species which are often in contact leading to 
a larger fraction of active links, see Fig. |3] 
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The central feature of the model is the ability to sustain 
high diversity, even in the limit a 0. Fig. [4] quantifies 
this for a = 200 x 200 system, which shows the finite 
diversity D in the limit of a — >■ for 7 < 0.055. In 
this figure we also show D{'j) obtained from a quasistatic 
simulation where new species are only introduced when 
there is no ongoing population dynamics. For 7 > 7^ « 
0.055 the quasistatic diversity is -0(7) = 1. 

The quasistatic version of the model reflects the bio- 
logically interesting limit where evolution rates are much 
slower than any dynamics associated to the species pop- 
ulations. The quasistatic simulation has two metastable 
states when 7 < 7c, one with a high but finite diversity 
and one absorbing state where D = 1. In the simula- 
tions, the state of high diversity is reached by starting 
from a state with large diversity. In simulating the qua- 
sistatic dynamics there will occasionally be long periods 
where several species compete dynamically for the same 
area. To shorten these periods we eliminate all outgoing 
links from one randomly chosen active species. If this 
does not stop the dynamics, the elimination procedure is 
repeated until the system eventually freezes. After the 
system has frozen, the eliminated links are re-introduced 
and the quasistatic simulation is continued by introduc- 
ing a new random species. The obtained steady state 
diversity (black line in Fig. |4]) is close to the full simula- 
tion at a = 0.01. 

Importantly, the obtained large diversity at 7 < 7c 
is also found in an extended model where species some- 
times compete for free space created by occasional death 
of individuals, i.e. the system behaviour is robust to ex- 
ternal pertubations of the type considered in |26) . In the 
quasistatic limit the sharp transition to high diversity is 
maintained even when 10% of the systems sites are killed 
prior to the introduction of each new species. For finite 
a, occasional death will soften the transition. 

In the insert of Fig. |4j we examine a variant of the qua- 
sistatic model where two new species are introduced si- 
multaneously instead of only one. Thereby the absorbing 
state at D — 1 disappears and the previously metastable 
state at high diversity becomes a true steady state. One 
sees that the obtained steady state diversity for this vari- 
ant of the model is close to that obtained from the a — 
limit of the standard model. Also, the insert in Fig. |4] 
shows that the behavior of the one- and two-species intro- 
duction models have a similar high diversity state when 
iV — 00. Both models exhibit a phase transition to the 
high diversity state for 7 < 7c w 0.055. The transition is 
maintained also when we allow the invasion rates to vary 
between the species, as well as when we implement the 
model on a triangular lattice. Close to the critical point, 
the finite but small a simulations exhibits pronounced 
bi-stability with rare transitions between a high and low 
diversity state. Transitions from low to high diversity are 
created through a state with many disconnected patches 
of a few species, whereas the breakdown of high diversity 
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FIG. 4. Diversity D as function of 7 for different a values. 
The sliaded areas mark the difference between the spatial 
model (upper bounding lines) and the random neighbor ver- 
sion (dashed lower lines). For the non-spatial version the 
diversity D « 11 x a°-^^ + 1 for L = 200 is independent of 
7. The black line shows the diversity achieved from the qua- 
sistatic limit. The inset shows simulations in the quasistatic 
limit where two species are simultaneously introduced when- 
ever a conhguration freezes. The black line is identical to 
the quasistatic limit shown in the main panel. The inset also 
shows that D /Area acts as an order parameter for the system. 

occurs when each species is only represented by a few 
patches. 

For real sessile ecosystems, 7 can be estimated from 
fraction of active boundaries, which have been reported 
down to a value of 2.5% for sponge-corals competition 
[31j . The low 7 is also consistent with observations on 
crustose lichens |30l . 

Discussion.- The model results in an interesting in- 
terplay between an interaction network |32| I33j and spa- 
tial dynamics. That is, the spatial configuration of the 
species limits the number of active links in the net- 
work and as a results most links are passive and have 
no influence on the dynamics. This is in contrast to 
the common network assumption that all links repre- 
sent ongoing interactions. In biological systems most 
links are transient. To explore further the influence of 
the space-network coupling, we have compared our spa- 
tial model with an equivalent random-neighbor version 
where lattice sites interact at random and independent 
of their spatial locations. The non-spatial model is a 
noisy version of average population equations of the form 
TT = r(z, - Tij,i)Xaj, with the addi- 

tional rule that species with near zero population are 
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removed. By comparing the dashed and sohd lines in 
Fig. |4j we see that the behavior of the random-neighbor 
model differs from our model, (i) D is independent of 7 
and in general much lower than for the spatial model. In 
the quasistatic limit then D = 1 for all 7 values, (ii) At 
moderate a and 7 the random neighbor model can de- 
velop sustained oscillatory states, resembling the hyper- 
cycles of Eigen and Shuster [1]. The size and duration 
of these cycles decrease with the rate of introduction of 
new species a. In the spatial model these cycles emerge in 
restricted spatial regions, cause redistribution of species 
boundaries and give rise to sub-populations separated by 
neutral species. 

The process that increases diversity in the model 
speaks to an interesting interplay between a sympatric 
and an allopatric speciation process. Although we have 
not implemented genealogical species relationships, the 
algorithm might be seen as a minimal evolution model, 
where new species are related to the species in the re- 
gion where they originate. In this framework, the spatial 
division of an existing species population by a neutral 
species and the subsequent invasion of one fragment by 
a new species then amount to an allopatric speciation 
event. Overall, our model can create two species in a 
region previously covered by one species. 

The main feature of our spatial model is the emergence 
of a robust ecological system with multiple co-existing 
species. The species diversity depends crucially on both 
the limitations imposed by the interaction network, and 
on the spatial positioning of the individual species. If 
either space or network constraint is absent the ecosys- 
tem complexity cannot be maintained. The co-existence 
is dynamic in the sense that it includes both extinction 
as well as occasional re-emergence of new species. In 
pS] it was found that non-transitive relationships in- 
deed prolonged co-existence of many species. However 
in that model a substantial fraction of the lattice was 
often cleared by external disturbances, and any of the 
initially present 10 species were repeatedly introduced on 
empty lattice positions. In contrast by allowing introduc- 
tion of species with new random interactions, our model 
demonstrated a self-organization towards an ecosystem 
with substantial diversity provided that the interaction 
probability 7 < 7c ^ 0.055. A system diversity which 
in fact remains even when we artificially cut cyclic re- 
lationships in the quasistatic limit, demonstrating that 
steady state species diversity crucially depends on chain- 
like spatial interactions that involves at least 3 species. 

In an evolutionary perspective, an interesting aspect of 
our model is the fact that it does not rely on any ultimate 
fitness landscape. At any point in space and time, the 
" fitness" of a species is dictated by its neighborhood [Ml 
155] . Our model supplements these studies by proposing 
a self-organized allopatric speciation mechanism which 
suggest a minimal sustainable diversity for ecosystems of 



mutually exclusive species in 2-dimensions. 
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